/********************************************************************************
 *    Copyright (C) 2019 Quark Matter Research Center, IMP,CAS			*
 *                                                                              *
 *              This software is distributed under the terms of the             *
 *              GNU Lesser General Public Licence (LGPL) version 3,             *
 *                  copied verbatim in the file "LICENSE"                       *
 ********************************************************************************/
// in root all sizes are given in cm




#include "TSystem.h"
#include "TGeoManager.h"
#include "TGeoVolume.h"
#include "TGeoMaterial.h"
#include "TGeoMedium.h"
#include "TGeoPgon.h"
#include "TGeoXtru.h"
#include "TGeoMatrix.h"
#include "TGeoCompositeShape.h"
#include "TFile.h"
#include "TString.h"
#include "TList.h"
#include "TROOT.h"



#include <iostream>

using namespace TMath;
const double cm=1.;
const double mm=0.1*cm;
const double um=0.001*mm;
const double pi=TMath::Pi();
const double deg=pi/180;

const double MWDC_thickness=110*mm; //thickness of MWDC.
const double MWDC_Al_cover_thickness=5*mm;
const int MWDC_Nwires_of_single_anode_layer= 80;               // number of anode wires in the single layer of MWDC, all are 80 for X, Xprim, V, Vprim, U, and Uprim layers
const double D_neighbouring_anode_wires = 10*mm;
const double D_cathode_and_anode_wires = 5*mm;
double MWDC_X_first_anode_wires_pos= 13*mm;
double MWDC_Xprim_first_anode_wires_pos= 17*mm;
double MWDC_V_first_anode_wires_pos= 124.8*mm;
double MWDC_Vprim_first_anode_wires_pos= 129.4*mm;
double MWDC_U_first_anode_wires_pos= 129.4*mm;
double MWDC_Uprim_first_anode_wires_pos=124.8*mm;


const double MWDC_outer_height=860*mm;
const double MWDC_outer_width=980*mm;
const double MWDC_inner_height=600*mm;
const double MWDC_inner_width=820*mm;

const double MWDC1_posX=60*cm;
const double MWDC1_posY=0.;
const double MWDC1_posZ=250*cm;
const double MWDC2_posX=100*cm;
const double MWDC2_posY=0.;
const double MWDC2_posZ=380*cm;
const double MWDC3_posX=140*cm;
const double MWDC3_posY=0.;
const double MWDC3_posZ=520*cm;
double MWDC1_FrontFace_posZ = MWDC1_posZ-MWDC_thickness/2.;
double MWDC2_FrontFace_posZ = MWDC2_posZ-MWDC_thickness/2.;
double MWDC3_FrontFace_posZ = MWDC3_posZ-MWDC_thickness/2.;

double Mylar_thickness = 1*mm;

const double thickness_of_single_board_layer=5*mm;
const double sensitive_layer_thickness= 2*thickness_of_single_board_layer;
////   
double LayerX_posZ_inMWDC= -8*thickness_of_single_board_layer;
double LayerXprim_posZ_inMWDC= -6*thickness_of_single_board_layer;
double LayerV_posZ_inMWDC= -1*thickness_of_single_board_layer;
double LayerVprim_posZ_inMWDC=1*thickness_of_single_board_layer;
double LayerU_posZ_inMWDC= 6*thickness_of_single_board_layer;
double LayerUprim_posZ_inMWDC= 8*thickness_of_single_board_layer;




// Name of geometry version and output file
TString geoPath = gSystem->Getenv("VMCWORKDIR");
const TString geoVersion = "mwdc_v1";
const TString FileName = geoPath+"/geometry/"+geoVersion + ".root";
const TString FileName1 = geoPath+"/geometry/"+geoVersion + "_geomanager.root";

// Names of the different used materials which are used to build the modules
// The materials are defined in the global media.geo file
const TString MWDCVolumeMedium  = "MWDCmixture";
const TString MWDC_mylar_str    = "mylar";
const TString Aluminum_str      = "Aluminum";
const TString PCB_str           = "pcbmvd";
const TString MWDCAnodeMedium   = "MWDC_anode";
const TString MWDCCathodeMedium = "MWDC_cathode";


// some global variables
TGeoManager* gGeoMan = NULL;  // Pointer to TGeoManager instance

// Forward declarations
void create_materials_from_media_file();
TLine calculate_wire_point(int N, TString plane, double x0, double w,double h);

void Create_MWDC_v1(){
	
	create_materials_from_media_file();
	
	gGeoMan=(TGeoManager*)gROOT->FindObject("FAIRGeom");
	gGeoMan->SetVisLevel(7);

	TGeoVolume* topMWDC= new TGeoVolumeAssembly("TOPMWDC");
	gGeoMan->SetTopVolume(topMWDC);

	/****** MWDC contains MWDC1, MWDC2 and MWDC3     */
	TGeoVolume* MWDC_array= new TGeoVolumeAssembly(geoVersion);
	topMWDC->AddNode(MWDC_array, 1);
	TGeoMedium* MWDCVolMed   = gGeoMan->GetMedium(MWDCVolumeMedium);
	TGeoMedium* MWDCAnodeMed   = gGeoMan->GetMedium(MWDCAnodeMedium);
	TGeoMedium* MWDCCathodeMed   = gGeoMan->GetMedium(MWDCCathodeMedium);
	TGeoMedium* MWDC_cover_med_mylar = gGeoMan->GetMedium(MWDC_mylar_str);  //前后foil的材料
	TGeoMedium* MWDC_cover_med_Al    = gGeoMan->GetMedium(Aluminum_str);
	TGeoMedium* Frame_med            = gGeoMan->GetMedium(PCB_str);

//	TGeoMedium* MWDCVolMed   = gGeoMan->GetMedium("air");
	TGeoRotation *zeroRotation= new TGeoRotation();
	zeroRotation->RotateX(0.);
	zeroRotation->RotateY(0.);
	zeroRotation->RotateZ(0.);
	TGeoCombiTrans* zero_CombiTrans=new TGeoCombiTrans(0,0,0,zeroRotation);
	
	TGeoRotation *gRot= new TGeoRotation();
	gRot->RotateX(0.);
	gRot->RotateY(15.);
	gRot->RotateZ(0.);
	/****** any MWDC has the same geometry ******/
	TGeoVolume* aMWDC= new TGeoVolumeAssembly("aMWDC");
	TGeoCombiTrans* mwdc1_combtrans=new TGeoCombiTrans(MWDC1_posX, MWDC1_posY, MWDC1_posZ,gRot);
	TGeoCombiTrans* mwdc2_combtrans=new TGeoCombiTrans(MWDC2_posX, MWDC2_posY, MWDC2_posZ,gRot);
	TGeoCombiTrans* mwdc3_combtrans=new TGeoCombiTrans(MWDC3_posX, MWDC3_posY, MWDC3_posZ,gRot);
	MWDC_array->AddNode(aMWDC, 1, mwdc1_combtrans);
	MWDC_array->AddNode(aMWDC, 2, mwdc2_combtrans);
	MWDC_array->AddNode(aMWDC, 3, mwdc3_combtrans);

	//double MWDC_gas_thickness=MWDC_thickness-2*MWDC_Al_cover_thickness;  //thickness of work gas 
	double MWDC_gas_thickness=MWDC_thickness;  //thickness of work gas 
	double MWDC_inner_thickness=MWDC_gas_thickness;    // thickness of PCB frame
	TGeoBBox* MWDC_Al_box=new TGeoBBox("MWDC_Al_box", MWDC_outer_width/2, MWDC_outer_height/2, MWDC_Al_cover_thickness/2);
	TGeoBBox* MWDC_Al_hole=new TGeoBBox("MWDC_Al_hole", MWDC_inner_width/2, MWDC_inner_height/2, MWDC_Al_cover_thickness/2+2);
	TGeoCompositeShape *MWDC_Al_shape= new TGeoCompositeShape("MWDC_Al_shape","MWDC_Al_box-MWDC_Al_hole");
	TGeoBBox* MWDC_PCB_box=new TGeoBBox("MWDC_PCB_box", MWDC_outer_width/2, MWDC_outer_height/2, MWDC_inner_thickness/2);
	TGeoBBox* MWDC_PCB_hole=new TGeoBBox("MWDC_PCB_hole", MWDC_inner_width/2, MWDC_inner_height/2, MWDC_inner_thickness/2+2);
	TGeoCompositeShape *MWDC_PCB_frame_shape= new TGeoCompositeShape("MWDC_PCB_frame_shape","MWDC_PCB_box-MWDC_PCB_hole");
	TGeoBBox* MWDC_gas_box=new TGeoBBox("MWDC_gas_box", MWDC_inner_width/2, MWDC_inner_height/2, MWDC_inner_thickness/2);

	TGeoVolume *MWDC_PCB_frame= new TGeoVolume("MWDC_PCB_frame", MWDC_PCB_frame_shape, Frame_med);
	MWDC_PCB_frame->SetLineColor(kYellow);
	MWDC_PCB_frame->SetTransparency(30);
	aMWDC->AddNode(MWDC_PCB_frame,1, zero_CombiTrans);
	
	TGeoVolume *MWDC_Al_cover = new TGeoVolume("MWDC_Al_cover", MWDC_Al_shape, MWDC_cover_med_Al);
	MWDC_Al_cover->SetLineColor(kGray);
	MWDC_Al_cover->SetTransparency(30);
	double Z_Al1= (MWDC_thickness-MWDC_Al_cover_thickness)/2;
	cout<<"Z_Al_cover: "<<Z_Al1/mm<<" mm"<<endl;
	TGeoCombiTrans* Al_CombiTrans1=new TGeoCombiTrans(0,0,-Z_Al1,zeroRotation);
	TGeoCombiTrans* Al_CombiTrans2=new TGeoCombiTrans(0,0,Z_Al1,zeroRotation); 
	aMWDC->AddNode(MWDC_Al_cover,1, Al_CombiTrans1);    // front Al cover
	aMWDC->AddNode(MWDC_Al_cover,2, Al_CombiTrans2);    //back Al cover

	TGeoVolume *MWDC_gas= new TGeoVolume("MWDC_gas", MWDC_gas_box, MWDCVolMed);
	MWDC_gas->SetLineColor(kGreen);
	//MWDC_gas->SetTransparency(10);
	aMWDC->AddNode(MWDC_gas,1, zero_CombiTrans);

	//// sensitive volumes in a MWDC, name as X, Xprim, V, Vprim, U, Uprim  ////
	//TGeoCombiTrans* X_CombiTrans=new TGeoCombiTrans(0,0,LayerX_posZ_inMWDC,zeroRotation);
	//TGeoCombiTrans* Xprim_CombiTrans=new TGeoCombiTrans(0,0,LayerXprim_posZ_inMWDC,zeroRotation);
	//TGeoCombiTrans* V_CombiTrans=new TGeoCombiTrans(0,0,LayerV_posZ_inMWDC,zeroRotation);
	//TGeoCombiTrans* Vprim_CombiTrans=new TGeoCombiTrans(0,0,LayerVprim_posZ_inMWDC,zeroRotation);
	//TGeoCombiTrans* U_CombiTrans=new TGeoCombiTrans(0,0,LayerU_posZ_inMWDC,zeroRotation);
	//TGeoCombiTrans* Uprim_CombiTrans=new TGeoCombiTrans(0,0,LayerUprim_posZ_inMWDC,zeroRotation);
  
	TGeoVolume *MWDC_sensitive_layerX= new TGeoVolumeAssembly("MWDC_sensitive_layerX");
	MWDC_gas->AddNode(MWDC_sensitive_layerX,1, new TGeoTranslation(0, 0, LayerX_posZ_inMWDC));
	cout<<"LayerX_posZ: "<<LayerX_posZ_inMWDC/mm<<" mm"<<endl;
	
	TGeoVolume *MWDC_sensitive_layerXprim= new TGeoVolumeAssembly("MWDC_sensitive_layerXprim");
	TGeoVolume *MWDC_sensitive_layerV= new TGeoVolumeAssembly("MWDC_sensitive_layerV");
	TGeoVolume *MWDC_sensitive_layerVprim= new TGeoVolumeAssembly("MWDC_sensitive_layerVprim");
	TGeoVolume *MWDC_sensitive_layerU= new TGeoVolumeAssembly("MWDC_sensitive_layerU");
	TGeoVolume *MWDC_sensitive_layerUprim= new TGeoVolumeAssembly("MWDC_sensitive_layerUprim");
	MWDC_gas->AddNode(MWDC_sensitive_layerXprim,1, new TGeoTranslation(0, 0, LayerXprim_posZ_inMWDC));
	MWDC_gas->AddNode(MWDC_sensitive_layerV,1, new TGeoTranslation(0, 0, LayerV_posZ_inMWDC));
	MWDC_gas->AddNode(MWDC_sensitive_layerVprim,1, new TGeoTranslation(0, 0, LayerVprim_posZ_inMWDC));
	MWDC_gas->AddNode(MWDC_sensitive_layerU,1, new TGeoTranslation(0, 0, LayerU_posZ_inMWDC));
	MWDC_gas->AddNode(MWDC_sensitive_layerUprim,1, new TGeoTranslation(0, 0, LayerUprim_posZ_inMWDC));

	////// wire module for X sensitive layer
	for(int WireID = 0;WireID<MWDC_Nwires_of_single_anode_layer;WireID++){
	TGeoXtru* WireUnit_shape = new TGeoXtru(2);    /// number of Z sections: 2
	double xx[4], yy[4];
	TLine line1=calculate_wire_point(WireID, "X",MWDC_X_first_anode_wires_pos-D_cathode_and_anode_wires,MWDC_inner_height, MWDC_inner_width);
	xx[0]= line1.GetX1();
	yy[0]= line1.GetY1();
	xx[3]= line1.GetX2();
	yy[3]= line1.GetY2();
		TLine line2=calculate_wire_point(WireID+1, "X",MWDC_X_first_anode_wires_pos-D_cathode_and_anode_wires,MWDC_inner_height, MWDC_inner_width);
	xx[1]= line2.GetX1();
	yy[1]= line2.GetY1();
	xx[2]= line2.GetX2();
	yy[2]= line2.GetY2();
	WireUnit_shape->DefinePolygon(4,xx,yy);
	
	// DefineSection(0, z0, x0, y0, scale0);   // Z position, offset and scale for first section
	// DefineSection(nz-1, zn, xn, yn, scalen); // Z position, offset and scale for first section
	WireUnit_shape->DefineSection(0, - thickness_of_single_board_layer, 0, 0, 1.0);
	WireUnit_shape->DefineSection(1, + thickness_of_single_board_layer, 0, 0, 1.0);
	TGeoVolume *WireUnit = new TGeoVolume(TString("MWDC_WireUnit_X")+Form("_%d",WireID), WireUnit_shape, MWDCVolMed);
	WireUnit->SetLineColor(kSpring);
	WireUnit->SetTransparency(90);
	MWDC_sensitive_layerX->AddNode(WireUnit,1,new TGeoTranslation(0, 0, 0));
	}

	////// wire module for Xprim sensitive layer
	for(int WireID = 0;WireID<MWDC_Nwires_of_single_anode_layer;WireID++){
	TGeoXtru* WireUnit_shape = new TGeoXtru(2);    /// number of Z sections: 2
	double xx[4], yy[4];
	TLine line1=calculate_wire_point(WireID, "X",MWDC_Xprim_first_anode_wires_pos-D_cathode_and_anode_wires,MWDC_inner_height, MWDC_inner_width);
	xx[0]= line1.GetX1();
	yy[0]= line1.GetY1();
	xx[3]= line1.GetX2();
	yy[3]= line1.GetY2();
		TLine line2=calculate_wire_point(WireID+1, "X",MWDC_Xprim_first_anode_wires_pos-D_cathode_and_anode_wires,MWDC_inner_height, MWDC_inner_width);
	xx[1]= line2.GetX1();
	yy[1]= line2.GetY1();
	xx[2]= line2.GetX2();
	yy[2]= line2.GetY2();
	
	WireUnit_shape->DefinePolygon(4,xx,yy);
	// DefineSection(0, z0, x0, y0, scale0);   // Z position, offset and scale for first section
	// DefineSection(nz-1, zn, xn, yn, scalen); // Z position, offset and scale for first section
	WireUnit_shape->DefineSection(0, - thickness_of_single_board_layer, 0, 0, 1.0);
	WireUnit_shape->DefineSection(1, + thickness_of_single_board_layer, 0, 0, 1.0);
	TGeoVolume *WireUnit = new TGeoVolume(TString("MWDC_WireUnit_Xprim")+Form("_%d",WireID), WireUnit_shape, MWDCVolMed);
	WireUnit->SetLineColor(kSpring);
	WireUnit->SetTransparency(90);
	MWDC_sensitive_layerXprim->AddNode(WireUnit,1,new TGeoTranslation(0, 0, 0));
	}

	////// wire module for V sensitive layer
	for(int WireID = 0;WireID<MWDC_Nwires_of_single_anode_layer;WireID++){
	TGeoXtru* WireUnit_shape = new TGeoXtru(2);    /// number of Z sections: 2
	double xx[4], yy[4];
	TLine line1=calculate_wire_point(WireID, "V",MWDC_V_first_anode_wires_pos-D_cathode_and_anode_wires,MWDC_inner_height, MWDC_inner_width);
	xx[0]= line1.GetX1();
	yy[0]= line1.GetY1();
	xx[3]= line1.GetX2();
	yy[3]= line1.GetY2();
		TLine line2=calculate_wire_point(WireID+1, "V",MWDC_V_first_anode_wires_pos-D_cathode_and_anode_wires,MWDC_inner_height, MWDC_inner_width);
	xx[1]= line2.GetX1();
	yy[1]= line2.GetY1();
	xx[2]= line2.GetX2();
	yy[2]= line2.GetY2();
	
	WireUnit_shape->DefinePolygon(4,xx,yy);
	// DefineSection(0, z0, x0, y0, scale0);   // Z position, offset and scale for first section
	// DefineSection(nz-1, zn, xn, yn, scalen); // Z position, offset and scale for first section
	WireUnit_shape->DefineSection(0, - thickness_of_single_board_layer, 0, 0, 1.0);
	WireUnit_shape->DefineSection(1, + thickness_of_single_board_layer, 0, 0, 1.0);
	TGeoVolume *WireUnit = new TGeoVolume(TString("MWDC_WireUnit_V")+Form("_%d",WireID), WireUnit_shape, MWDCVolMed);
	WireUnit->SetLineColor(kSpring);
	WireUnit->SetTransparency(90);
	MWDC_sensitive_layerV->AddNode(WireUnit,1,new TGeoTranslation(0, 0, 0));
	}

	////// wire module for Vprim sensitive layer
	for(int WireID = 0;WireID<MWDC_Nwires_of_single_anode_layer;WireID++){
	TGeoXtru* WireUnit_shape = new TGeoXtru(2);    /// number of Z sections: 2
	double xx[4], yy[4];
	TLine line1=calculate_wire_point(WireID, "V",MWDC_Vprim_first_anode_wires_pos-D_cathode_and_anode_wires,MWDC_inner_height, MWDC_inner_width);
	xx[0]= line1.GetX1();
	yy[0]= line1.GetY1();
	xx[3]= line1.GetX2();
	yy[3]= line1.GetY2();
		TLine line2=calculate_wire_point(WireID+1, "V",MWDC_Vprim_first_anode_wires_pos-D_cathode_and_anode_wires,MWDC_inner_height, MWDC_inner_width);
	xx[1]= line2.GetX1();
	yy[1]= line2.GetY1();
	xx[2]= line2.GetX2();
	yy[2]= line2.GetY2();
	
	WireUnit_shape->DefinePolygon(4,xx,yy);
	// DefineSection(0, z0, x0, y0, scale0);   // Z position, offset and scale for first section
	// DefineSection(nz-1, zn, xn, yn, scalen); // Z position, offset and scale for first section
	WireUnit_shape->DefineSection(0, - thickness_of_single_board_layer, 0, 0, 1.0);
	WireUnit_shape->DefineSection(1, + thickness_of_single_board_layer, 0, 0, 1.0);
	TGeoVolume *WireUnit = new TGeoVolume(TString("MWDC_WireUnit_Vprim")+Form("_%d",WireID), WireUnit_shape, MWDCVolMed);
	WireUnit->SetLineColor(kSpring);
	WireUnit->SetTransparency(90);
	MWDC_sensitive_layerVprim->AddNode(WireUnit,1,new TGeoTranslation(0, 0, 0));
	}

	////// wire module for U sensitive layer
	for(int WireID = 0;WireID<MWDC_Nwires_of_single_anode_layer;WireID++){
	TGeoXtru* WireUnit_shape = new TGeoXtru(2);    /// number of Z sections: 2
	double xx[4], yy[4];
	TLine line1=calculate_wire_point(WireID, "U",MWDC_U_first_anode_wires_pos-D_cathode_and_anode_wires,MWDC_inner_height, MWDC_inner_width);
	xx[0]= line1.GetX1();
	yy[0]= line1.GetY1();
	xx[3]= line1.GetX2();
	yy[3]= line1.GetY2();
		TLine line2=calculate_wire_point(WireID+1, "U",MWDC_U_first_anode_wires_pos-D_cathode_and_anode_wires,MWDC_inner_height, MWDC_inner_width);
	xx[1]= line2.GetX1();
	yy[1]= line2.GetY1();
	xx[2]= line2.GetX2();
	yy[2]= line2.GetY2();
	
	WireUnit_shape->DefinePolygon(4,xx,yy);
	// DefineSection(0, z0, x0, y0, scale0);   // Z position, offset and scale for first section
	// DefineSection(nz-1, zn, xn, yn, scalen); // Z position, offset and scale for first section
	WireUnit_shape->DefineSection(0, - thickness_of_single_board_layer, 0, 0, 1.0);
	WireUnit_shape->DefineSection(1, + thickness_of_single_board_layer, 0, 0, 1.0);
	TGeoVolume *WireUnit = new TGeoVolume(TString("MWDC_WireUnit_U")+Form("_%d",WireID), WireUnit_shape, MWDCVolMed);
	WireUnit->SetLineColor(kSpring);
	WireUnit->SetTransparency(90);
	MWDC_sensitive_layerU->AddNode(WireUnit,1,new TGeoTranslation(0, 0, 0));
	}

	////// wire module for Uprim sensitive layer
	for(int WireID = 0;WireID<MWDC_Nwires_of_single_anode_layer;WireID++){
	TGeoXtru* WireUnit_shape = new TGeoXtru(2);    /// number of Z sections: 2
	double xx[4], yy[4];
	TLine line1=calculate_wire_point(WireID, "U",MWDC_Uprim_first_anode_wires_pos-D_cathode_and_anode_wires,MWDC_inner_height, MWDC_inner_width);
	xx[0]= line1.GetX1();
	yy[0]= line1.GetY1();
	xx[3]= line1.GetX2();
	yy[3]= line1.GetY2();
		TLine line2=calculate_wire_point(WireID+1, "U",MWDC_Uprim_first_anode_wires_pos-D_cathode_and_anode_wires,MWDC_inner_height, MWDC_inner_width);
	xx[1]= line2.GetX1();
	yy[1]= line2.GetY1();
	xx[2]= line2.GetX2();
	yy[2]= line2.GetY2();
	
	WireUnit_shape->DefinePolygon(4,xx,yy);
	// DefineSection(0, z0, x0, y0, scale0);   // Z position, offset and scale for first section
	// DefineSection(nz-1, zn, xn, yn, scalen); // Z position, offset and scale for first section
	WireUnit_shape->DefineSection(0, - thickness_of_single_board_layer, 0, 0, 1.0);
	WireUnit_shape->DefineSection(1, + thickness_of_single_board_layer, 0, 0, 1.0);
	TGeoVolume *WireUnit = new TGeoVolume(TString("MWDC_WireUnit_Uprim")+Form("_%d",WireID), WireUnit_shape, MWDCVolMed);
	WireUnit->SetLineColor(kSpring);
	WireUnit->SetTransparency(90);
	MWDC_sensitive_layerUprim->AddNode(WireUnit,1,new TGeoTranslation(0, 0, 0));
	}

	cout<<"Voxelizing."<<endl;
	topMWDC->Voxelize("");
	gGeoMan->CloseGeometry();
	topMWDC->Draw();
	gGeoMan->CheckOverlaps(0.001);
	gGeoMan->PrintOverlaps();
	gGeoMan->Test();

	TFile* outfile = TFile::Open(FileName,"RECREATE");
	topMWDC->Write();
	outfile->Close();

	TFile* outfile1 = TFile::Open(FileName1,"RECREATE");
	gGeoMan->Write();
	outfile1->Close();
}

void create_materials_from_media_file()
{
  // Use the FairRoot geometry interface to load the media which are already defined
  FairGeoLoader* geoLoad = new FairGeoLoader("TGeo", "FairGeoLoader");
  FairGeoInterface* geoFace = geoLoad->getGeoInterface();
  //TString geoPath = gSystem->Getenv("VMCWORKDIR");
  TString geoFile = geoPath + "/geometry/media.geo";
  geoFace->setMediaFile(geoFile);
  geoFace->readMedia();
  // Read the required media and create them in the GeoManager
  FairGeoMedia* geoMedia = geoFace->getMedia();
  FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder();
  FairGeoMedium* air              = geoMedia->getMedium("air");
  FairGeoMedium* MWDCmixture     = geoMedia->getMedium("MWDCmixture");
  FairGeoMedium* MWDC_Anode     = geoMedia->getMedium("MWDC_anode");
  FairGeoMedium* MWDC_Cathode     = geoMedia->getMedium("MWDC_cathode");
	FairGeoMedium* MWDC_mylar     = geoMedia->getMedium("mylar");
	FairGeoMedium* Aluminum     = geoMedia->getMedium("Aluminum");
	FairGeoMedium* PCB = geoMedia->getMedium("pcbmvd");

  // include check if all media are found

  geoBuild->createMedium(air);
  geoBuild->createMedium(MWDCmixture);
  geoBuild->createMedium(MWDC_Anode);
  geoBuild->createMedium(MWDC_Cathode);
  geoBuild->createMedium(MWDC_mylar);
  geoBuild->createMedium(Aluminum);
  geoBuild->createMedium(PCB);
}

//void calculate_wire_point(int N, TString plane, double x0, double w,double h, double &x1,double &y1,double &x2,double &y2)
TLine calculate_wire_point(int N, TString plane, double x0, double h,double w)
{
/////////////////////////////////////////////////////////////////////////////////////////
/// calculate the points for the N wire crossing with the inner edges of the PCB frame //
//  N: the wireID, range from 0 to 80                                                  //  
//  plane: "X", "U", "V"                                                               //  
//  x0: the distance from the first wire to the nearest inner pcb edge.                //  
//  w:  width of the inner frame                                                       //  
//  h:  height of the inner frame                                                      //  
//  (x1,y1), (x2,y2) are the up and down points for the wire crosing the inner frame   //  
/////////////////////////////////////////////////////////////////////////////////////////
	  // y = k(x-x0)+y0 , y1>y2  y0==h/2
		
	TLine line;
	if(N<0||N>80) {
		cout<<"N = "<<N<<" is not correct, N range should be from 0 to 80. "<<endl;
	return line;
	}
	double x1,y1,x2,y2;
	if(plane.EqualTo("X"))     // X plane
	{
		x1=-w/2+x0+N*D_neighbouring_anode_wires;
		y1=h/2;
		x2=x1;
		y2=-h/2;

	}
	if(plane.EqualTo("V"))     // V plane
	{
		double theta=30*deg;
		double Xn=-w/2+x0+N*D_neighbouring_anode_wires/cos(theta);
		double YL=h/2+1/TMath::Tan(theta)*(-w/2-Xn);
		double YR=h/2+1/TMath::Tan(theta)*(w/2-Xn);
		double XD=-h*TMath::Tan(theta)+Xn;

		x1=TMath::Min(Xn, w/2);
		y1=TMath::Min(h/2, YR);  
		x2=TMath::Max(-w/2, XD);
		y2=TMath::Max(-h/2, YL);

		//cout<<"cos(30deg): "<<cos(theta)<<", Xn: "<<Xn<<", XD: "<<XD<<",  YL: "<<YL<<", YR: "<<YR<<endl;

	}
	if(plane.EqualTo("U"))     // U plane
	{
		double theta=-30*deg;
		double Xn=-w/2+x0+N*D_neighbouring_anode_wires/cos(theta);
		double YL=-h/2+1/TMath::Tan(theta)*(-w/2-Xn);
		double YR=-h/2+1/TMath::Tan(theta)*(w/2-Xn);
		double XUp=h*TMath::Tan(theta)+Xn;

		x1=TMath::Max(-w/2, XUp);
		y1=TMath::Min(h/2, YL);  
		x2=TMath::Min(Xn, w/2);
		y2=TMath::Max(-h/2, YR);

		//cout<<"cos(30deg): "<<cos(theta)<<", Xn: "<<Xn<<", XD: "<<XD<<",  YL: "<<YL<<", YR: "<<YR<<endl;

	}
	line.SetX1(x1);
	line.SetY1(y1);
	line.SetX2(x2);
	line.SetY2(y2);
	//line.Print();
	return line;
}

